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ABSTRACT 


It  is  well  recognized  that  the  material  properties 
in  fiber  reinforced  components  are  strongly  dependent  on 
the  fiber  orientation.  In  mold  filling  processes  involving 
short  fiber  reinforced  composites,  fiber  orientation  occurs 
as  a  result  of  the  flow  induced  stresses.  It  is  important 
to  be  able  to  predict  this  flow  induced  orientation. 

A  numerical  method  has  been  developed  previously 
to  predict  the  in-plane  fiber  orientation  in  plane  flow. 
This  scheme  is  refined  to  enable  predictions  for  fiber 
orientation  in  axisymmetric  flow.  The  numerical  method 
is  verified  by  comparing  numerical  and  analytical  solu¬ 
tions  for  fiber  orientation  in  Poiseuille  flow. 

The  fiber  orientation  may  be  described  by  certain 
orientation  parameters,  which  relate  the  degree  of  col- 
limation  to  the  material  properties.  These  orientation 
parameters  are  incorporated  into  both  the  plane  and  axi¬ 
symmetric  flow  algorithms,  thus  providing  a  link  to 
available  structural  analysis  routines. 


lil 


Several  numerical  examples  of  practical 
importance  are  presented.  A  prediction  for  the  fiber 
orientation  near  molded  holes  is  established  by  determin¬ 
ing  the  fiber  orientation  resulting  from  flow  around  a 
circular  inclusion.  In  a  further  example,  the  process  for 
molding  an  axisymmetric  disk  is  simulated  to  determine  the 
numerically  predicted  fiber  orientation. 
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CHAPTER  1 .  INTRODUCTION 


Composite  materials  are  becoming  widely  used, 
especially  for  aerospace  applications  where  stiff,  strong 
lightweight  materials  are  required.  However,  traditional 
hand  lay-up  methods  for  fabricating  composite  materials 
are  expensive  and  time-consuming  and,  hence,  can  only  be 
used  when  increased  structural  performance  outweighs  cost 
requirements.  With  the  advent  of  less  expensive  fabrica¬ 
tion  techniques ,  composites  are  becoming  sound  economic 
alternatives  for  many  other  applications. 

One  promising  inexpensive  fabrication  method  is 
injection  molding,  a  process  in  which  a  plasticized 
charge  is  forced  under  pressure  into  the  cavity  of  a 
closed  die  where  it  (the  charge)  is  formed  into  the 
shape  of  that  cavity.  During  this  mold  filling  process, 
fiber  alignment  occurs  due  to  fluid  stresses  induced  by 
the  flow.  Since  the  material  properties  and  strength 
of  the  resulting  component  part  are  strongly  dependent 
on  the  orientation  of  the  fibers,  it  is  essential  to  be 
able  to  predict  this  flow  induced  fiber  orientation. 
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An  early  attempt  to  predict  flow-induced 
orientation  was  performed  by  Jeffery  [9].  Jeffery  devel¬ 
oped  equations  describing  the  motion  of  a  single  ellipsoi 
(Ja.1  particle  immersed  in  a  Newtonian  fluid.  Goldsmith  and 
Mason  [7]  have  used  Jeffery's  equations  to  solve  for  the 
orientation  of  ellipsoidal  particles  in  both  Couette  and 
hyperbolic  radial  flows.  Based  on  Goldsmith  and  Mason's 
predictions  for  orientation  in  hyperbolic  radial  flow, 
Goettler  et  al. [6]  presented  a  new  manufacturing  concept 
for  producing  short  fiber  reinforcement  in  extruded  rubber 
hoses  in  a  one-step  fabrication  process.  Their  idea  was  to 
alter  the  die  design  such  that  a  constriction  of  the  flow 
occurs  at  some  intermediate  point,  followed  by  an  expansion 
to  form  the  dimensions  of  the  final  product.  In  the  region 
of  expansion,  fiber  orientation  occurs  in  the  hoop  direc¬ 
tion,  thus  providing  the  necessary  reinforcement.  The  work 
by  Goettler  et  al.  represents  one  of  the  early  attempts  to 
control  the  orientation  in  short  fiber  molded  components 
based  on  analytical  techniques. 

Givler  [5]  has  developed  a  numerical  method  to 
predict  in-plane  fiber  orientation  in  plane  flow.  The 
method  consists  of  solving  the  flow  equations  via  the 
finite  element  method  and,  subsequently,  numerically  inte¬ 
grating  Jeffery's  equations  to  determine  fiber  orientation. 
Givler' s  work  represents  a  first  attempt  to  develop  a 
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general  method  for  determining  fiber  orientation,  although 
it  is  restricted  to  in-plane  fibers  in  plane  flow. 

Once  the  fiber  orientation  is  determined,  material 
properties  may  be  predicted.  Using  the  "aggregate  model," 
McCullough  et  al.  [13]  have  shown  that  two  orientation 
parameters  are  needed  to  relate  the  degree  of  orientation 
to  the  material  properties  for  planar  fiber  orientation 
distributions,  if  the  orientation  distribution  is  symmetric 
about  the  mode  of  the  distribution.  This  mode  orientation 
angle  serves  to  isolate  the  local  principal  material  axes. 

The  first  objective  of  this  work  is  to  determine 
the  orientation  parameters  and  mean  orientation  angle  for 
the  in-plane  fiber  distributions  in  plane  flow  from 
Givler's  numerical  orientation  scheme,  thus  linking  Givler's 
algorithm  to  available  structural  analysis  programs  to  form 
a  complete  numerical  system  capable  of  evaluating  the  in¬ 
fluence  of  mold  designs  on  the  performance  of  molded  parts. 

The  second  objective  of  this  work  is  to  determine 
fiber  orientation  in  axisymmetric  flow.  In  the  axisymmet- 
ric  development,  the  assumption  of  planar  orientation  is 
lifted.  Again,  the  ultimate  goal  is  to  portray  the  fiber 
orientation  in  terms  of  orientation  parameters.  A  new  set 
of  orientation  parameters  is  needed  in  view  of  the  non- 
planar  distribution  assumption.  McGee  [14]  has  determined 
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that  four  orientation  parameters  are  necessary  to  compute 
the  material  properties  for  non-planar  distributions.  It 
is  these  four  orientation  parameters  which  are  determined 
in  the  axisymmetric  scheme. 

The  contents  of  this  thesis  are  organized  in  the 
following  manner:  Chapter  2  introduces  Jeffery's  orienta¬ 
tion  equations  and  briefly  summarizes  the  various  assump¬ 
tions  made  in  the  development  of  these  equations.  The 
analysis  involved  in  determining  fiber  orientation  from 
Jeffery's  equations  is  then  explained.  Chapter  3  reviews 
Givler's  numerical  method  for  determining  the  in-plane 
fiber  orientation  in  plane  flow,  and  introduces  algorithms 
for  determining  the  mean  orientation  angle  and  the  orienta¬ 
tion  parameters.  Several  numerical  examples  are  presented. 
Chapter  4  develops  the  numerical  method  for  determining 
fiber  orientation  in  axisymmetric  flow  and,  subsequently, 
devises  the  scheme  for  computing  the  non-planar  orientation 
parameters.  As  in  Chapter  3,  several  numerical  examples 
are  included.  Results  and  conclusions,  as  well  as 
suggestions  for  further  research,  are  presented  in 
Chapter  5. 


CHAPTER  2.  GENERAL  ORIENTATION  EQUATIONS 

In  this  chapter,  orientation  equations  developed  by 
Jeffery  [9]  are  incorporated  to  solve  for  fiber  orientation. 
Jeffery  formulated  equations  to  determine  the  motion  of  a 
single  rigid  ellipsoidal  particle  immersed  in  a  viscous 
fluid  subject  to  the  following  assumptions: 

1.  Apart  from  the  local  disturbance  near  the  particle, 
the  fluid  motion  is  steady  and  varies  in  space  on  a 
scale  that  is  large  compared  with  the  dimensions  of 
the  particle. 

2.  The  fluid  which  surrounds  the  particle  is  incom¬ 
pressible  and  Newtonian. 

3.  The  fluid  velocity  is  low;  hence,  inertia  terms  may 
be  neglected  (creeping  flow) . 

4.  The  particle  is  non- sedimenting. 

A  summary  of  Jeffery's  analysis  is  presented  here. 

For  details,  the  interested  reader  is  referred  to  the 
original  reference. 

One  important  conclusion  drawn  by  Jeffery  is  that 
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under  the  above  assumptions,  the  particle  ultimately 
attains  the  velocity  of  the  fluid  which  it  displaces.  This 
enables  one  to  track  the  paths  taken  by  particles  simply  by 
determining  the  flow  streamlines. 


To  solve  for  the  particle  orientation,  a 

_  0  0  0 

rectangular  Cartesian  coordinate  set  of  axes  x^  ,  x2  ,  , 

fixed  in  the  particle,  is  defined  such  that  the  surface  of 
the  ellipsoid  is  described  by  the  expression: 


(x°) 2  (x°)2  _  (x°)2 
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It  is  natural  to  determine  the  orientation  of  these  axes 
relative  to  a  set  of  axes  x^  ,  x2  ,  x2  fixed  in  direction 
whose  origin  lies  at  the  center  of  the  particle.  The  rela¬ 
tive  orientation  between  the  two  sets  of  axes  is  described 
by  the  three  Euler  angles.  If  the  particle  is  an  ellipsoid 
of  revolution,  however,  only  two  of  the  three  Euler  angles 
are  needed  to  fully  describe  its  orientation,  as  shown  in 
Figure  2.1.  A  fiber  may  be  modelled  as  an  ellipsoid  of 
revolution. 


Jeffery  found  expressions  relating  the  hydrodynamic 
torque  acting  on  the  particle  to  the  spins  about  its  axes. 
With  the  assumption  of  creeping  flow,  the  torque  vanishes, 
and  by  relating  the  spins  about  the  axes  to  the  Euler  angles 
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It  is  desirable  to  determine  the  orientation 
relative  to  an  inertial  reference  frame  since  the  fluid 

kinematics  are  most  naturally  defined  relative  to  a  fixed 
reference  frame.  For  this  purpose  an  inertial  Cartesian 
coordinate  system  x^  ,  X£  ,  x^  is  introduced  such  that  each 
x.^  axis  points  in  the  same  direction  as  the  corresponding 
x.  axis.  Obviously,  the  orientation  of  the  particle  axes 
relative  to  the  inertial  frame  can  be  described  by  the  same 
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One  final  limitation  of  the  above  equations  is  that 
they  were  developed  for  a  single  particle;  hence,  they  do 
not  account  for  particle  interactions  in  a  fluid  containing 
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many  particles.  Indeed,  in  typical  mold  filling 
operations,  fiber  volume  fractions  of  30-40%  are  common 
and  fiber  interations  may  become  very  significant.  However, 
the  analysis  involved  in  dealing  with  such  interactions  is 
at  the  present  time  intractable?  therefore,  such  inter¬ 
actions  are  not  treated  in  this  work. 

2.1  Use  of  Jeffery's  Equations  to  Determine  Fiber 

Orientation 

Upon  inspection  of  (2.3)  and  (2.4),  one  notices 
the  existence  of  vorticity  and  rate  of  deformation  com¬ 
ponents.  These  components  must  first  be  determined  and 
substituted  into  Jeffery's  equations,  which  are  then 
integrated  to  obtain  the  fiber  orientation. 

The  solution  to  the  flow  equations  must  satisfy 
conservation  of  mass  and  momentum  requirements,  along  with 
the  constitutive  approximation  for  the  fluid.  Conservation 
of  mass  and  momentum  are  universal  laws  applicable  to  all 
fluids,  but  the  constitutive  relation  is  merely  an  approxi¬ 
mation  for  modelling  the  behavior  of  each  specific  fluid. 

In  this  work,  the  "fluid"  consists  of  a  Newtonian  fluid 
containing  a  concentration  of  suspended  fibers.  For  a  low 
volume  fraction  of  fibers,  the  influence  of  the  fibers  on 
the  behavior  of  the  overall  fluid  is  minimal ,  and  the 
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Newtonian  constitutive  assumption  is  valid.  However,  a 
larger  concentration  of  fibers  exerts  a  considerable  in¬ 
fluence  on  the  constitutive  behavior  of  the  suspension. 
Maschmeyer  and  Hill  [10]  have  shown  that  when  a  high  concen¬ 
tration  of  3  mm  fibers  is  mixed  in  a  Newtonian  fluid,  the 
suspension  becomes  highly  pseudoplastic.  Thus  a  power  law 
constitutive  assumption  appears  valid  for  suspensions  of 
high  concentration. 

The  degree  of  orientation  of  the  fibers  can  also 
be  expected  to  have  an  influence  on  the  constitutive  prop¬ 
erties  of  the  suspension.  One  certainly  anticipates  dif¬ 
ferent  constitutive  behavior  in  regions  where  fibers  are 
highly  aligned  than  in  regions  of  random  orientation.  To 
include  this  effect,  however,  the  relationship  between 
viscosity  and  degree  of  orientation  must  first  be  deter¬ 
mined.  Thereafter,  an  iterative  procedure  is  needed  since 
the  fluid  mechanics  has  an  effect  on  fiber  orientation 
which  in  turn  alters  the  fluid  mechanics.  However,  for 
the  remainder  of  this  work,  the  vicosity  dependence  on 
fiber  orientation  is  neglected. 


CHAPTER  3.  PLANE  FLOW 


This  chapter  deals  with  in-plane  fiber  orientation 
in  plane  flow.  Numerically,  the  flow  equations  are 
easier  to  solve  in  plane  flow  since  the  flow  is  two- 
dimensional.  Also,  the  fiber  orientation  equations  sim¬ 
plify  considerably,  and  in  special  cases,  can  be  solved 
analytically.  In  Section  3.1,  the  appropriate  simplifica¬ 
tions  are  introduced  into  the  fiber  orientation  equations. 
Furthermore,  it  is  shown  that  in-plane  fibers  maintain 
their  in-plane  orientation  and  one  of  the  orientation 
equations  is  eliminated.  Section  3.2  present  an  analyti¬ 
cal  solution  for  fiber  orientation  in  plane  Poiseuille 
flow.  This  solution  provides  a  basis  for  checking  the 
numerical  scheme.  The  numerical  method  is  developed  in 
Section  3.3.  The  orientation  parameters  which  relate  the 
degree  of  orientation  to  the  material  properties  are  in¬ 
troduced  in  Section  3 . 4  and  a  method  for  determining  them 
as  well  as  the  mode  of  the  distribution  in  terms  of  the 
orientations  of  a  finite  number  of  fibers  is  presented. 
Finally,  in  Section  3.5,  a  number  of  numerically  deter¬ 
mined  solutions  are  presented.  To  check  the  numerical 
scheme,  the  fiber  orientation  in  plane  Poiseuille  flow  is 
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evaluated  and  the  results  are  compared  with  the  analytical 
solution  in  Section  3.2.  Also,  the  fiber  orientation  in 
flow  around  a  circular  inclusion  in  both  finite  and  in¬ 
finite  width  channels  is  determined.  The  latter  examples 
are  of  pragmatic  importance  because  they  provide  predic¬ 
tions  for  the  orientation  around  molded  holes. 

3 . 1  Orientation  Equation  for  In-Plane  Fibers  in 

Plane  Flow 

Consider  plane  flow  in  which  variations  in  the 
x^  coordinate  direction  may  be  ignored.  For  this  type  of 
flow,  several  of  the  vorticity  and  rate  of  deformation 
tensor  components  are  identically  zero : 


0  0  0 


0  d 


22 


d 


23 


Inserting  the  appropriate  simplifications  into  Equations 
(2.3)  and  (2.4)  one  obtains: 
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D<J> 

jj£-=z1  +  Btd23  cos2<))1  -  Js(d22-d33)sin2cj)1]  (3.1) 

D01  B 

=  -^[2d23  sin2(f>^  sin203  +  (d22~d33)  cos2<{)^  sin203 

+  3 (d22+d33) sin203]  (3.2) 

Equation  (3.1)  clearly  indicates  that  the  in-plane  response 
of  a  fiber  is  independent  of  the  out-of-plane  orientation. 
If  one  considers  fibers  initially  oriented  in  the  plane 
of  the  flow  (0j=tT/2)  ,  then  equation  (3.2)  reduces  to  jjjr—  =  0 
which  indicates  that  a  fiber  remains  in-plane.  The  re¬ 
mainder  of  this  chapter  deals  exclusively  with  in-plane 
fibers;  hence,  only  the  single  differential  equation  (3.1) 
needs  to  be  solved  to  completely  determine  the  fiber 
orientation. 

It  is  worthwhile  to  note  for  later  reference  that, 
under  the  above  simplifications,  Equations  (2.1)  and  (2.2) 
reduce  to  the  single  differential  equation: 

_  _  _ 

=  z1  +  B[d23  cos2*1  -  h (d22“d33)sin2^1l  (3.3) 

Although  the  above  assumptions  have  considerably 
simplified  the  analysis,  (3.1)  and  (3.3)  are  still 


15 


extremely  difficult  to  solve  analytically  except  for  some 
very  simple  flows.  The  next  section  deals  with  a  flow 
where  an  analytical  solution  does  exist.  Thereafter,  a 
procedure  is  developed  for  numerically  determing  fiber 
orientation  for  more  complicated  flows. 

3.2  Fiber  Orientation  in  Plane  Poiseuille  Flow  - 

Analytical  Solution 

Consider  the  pressure-driven,  steady, creeping 
flow  of  an  incompressible  Newtonian  fluid  through  a  rec¬ 
tangular  channel  of  very  large  aspect  ratio  as  shown  in 
Figure  3. 1* where  the  flow  is  in  the  x3  direction.  The 
ratios  H/w  and  H/L  are  small  compared  to  unity  and  thus 
regions  away  from  the  edges  of  the  channel ,  the  veloc¬ 
ity  variations  occur  only  in  the  x3  direction.  The  fiber 

orientation  in  this  fully  developed  region  is  to  be 
determined . 

Under  conditions  of  unit  flow,  the  normalized 
velocity  profile  is 
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0  0  0 

o  o  3sr 

o  nr  o 


du_  2 

where  r  =  =  -3x2/H 

2 

Substituting  into  Equation  (3.1),  one  obtains: 


3^ 

3t  "  +  U33x3 


— ^ — (r  2cos2<j)1  +  sin2<j>1) 
rp  +1 


(3.4) 


If  one  deals  with  fibers  entering  the  domain  perpendicular 

3<f>i 

to  the  flow  streamline  (ie,  4>^=0.  at  x3=0) ,  then  ^  =  0  and 

(3.4)  reduces  to 


d<J>1  ri  9  2  2 

-5— -  =  - ^ -  (r  cos  +  sin  $.) 

to3  u3(rp2+1)  p  1  1 


and  the  resulting  solution  is 


rx3/u3 

tan  <j>^  =  r  tan  -  +1/r 

P  P 


or 
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tan  <j>,  =  r  tan 

J.  £r 

Assuming  a  fiber  aspect  ratio  of  50,  the  fiber  orientation 
is  plotted  in  Figure  3.2  for  H=l.  One  immediate  observa¬ 
tion  to  be  made  from  the  figure  is  that  a  distinct  layer 
of  fibers  aligned  parallel  to  the  streamlines  exists  near 
the  wall  boundary.  This  alignment  correlates  well  with  ex¬ 
perimental  observations  of  fiber  alignment  near  wall  bound¬ 
aries  in  short- fiber  injection  molded  components. 

Having  attained  an  analytical  solution,  one  may 
now  proceed  to  develop  a  numerical  scheme  to  predict  fiber 
orientation.  The  known  analytic  solution  in  this  section 
can  be  used  to  check  the  accuracy  of  the  numerical  method. 

3. 3  Numerical  Method 

Except  for  very  few  simple  flow  geometries. 

Equation  (3.1)  cannot  be  solved  by  analytical  methods. 
However,  one  may  wish  to  determine  fiber  orientation  for 
complicated  flow  situations  which  exist  in  many  molding 
processes.  With  this  motivation,  the  necessity  to  develop 
a  numerical  scheme  to  determine  fiber  orientation  is  ob¬ 
vious.  The  majority  of  the  material  covered  in  this  sec¬ 
tion  was  previously  developed  by  Givler  [5] . 


2x2x3 


(H2-x2 ) (rp+l/rp) 
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Figure  3.2  Fiber  orientation  in  plane  Poiseuille  flow  for  fibers  initially 
perpendicular  to  flow  streamlines 
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As  mentioned  previously,  the  first  step  in  solving 
for  fiber  orientation  is  to  obtain  the  flow  field.  . 

The  finite  element  method  has  been  used  extensively  in 
the  literature  to  solve  fluid  mechanics  problems  with  ex¬ 
cellent  results  and  consequently,  is  incorporated  in  this 
work.  The  details  of  the  method  employed  herein  are  fully 
described  in  References  2  and  5.  The  analysis . utilizes 
the  six-node  triangular  element  shown  in  Figure  3.3  to 
solve  the  plane  and,  later,  axisymmetric  flows.  In  the 
finite  element  analysis,  the  continuity,  momentum,  and 
constitutive  equations  are  reduced  to  a  set  of  algebraic 
equations  via  the  Galerkin  method  and  after  solving  this 
system  of  equations,  a  discretized  solution  is  obtained. 

From  the  discretized  solution  for  the  fluid 
mechanics,  it  remains  to  determine  the  fiber  orientation. 
This  is  accomplished  by  determining  the  orientation  of 
individual  fibers  as  they  traverse  the  flow  streamlines. 


3 


I  4  2 


Figure  3.3  Six- node  triangular  element 
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Obviously,  unless  some  collimation  is  evident,  one 
cannot  know  precisely  the  initial  orientations  of  all  fi~ 
bers.  Indeed,  if  the  initial  orientations  were  known,  it 
would  not  be  feasible  to  track  each  fiber  through  the  do¬ 
main.  An  alternative  approach  is  to  select  a  number  of 
streamlines  which  sensibly  cover  the  domain  in  which  the  fi¬ 
ber  orientation  is  to  be  determined.  Then  a  finite  number 
of  fibers  of  varying  initial  orientations  are  allowed  to 
traverse  each  streamline,  and  the  orientation  of  each  fiber 
is  computed  at  various  locations  along  the  streamline.  Thus 
information  is  obtained  about  the  dispersion  of  fibers  as 
well  as  the  principal  direction  of  orientation.  Usually, 
an  initially  random  orientation  distribution  is  input.  A 
random  distribution  will  be  defined  in  a  later  section. 

To  find  the  streamlines ,  one  must  first  determine 
the  stream  function,  ip .  From  the  numerical  solution, 
the  stream  function  values  are  known  at  the  nodal 
points  as  shown  in  Figure  3.4.  To  isolate  a  streamline, 
the  finite  element  grid  is  scanned  element  by 
element  to  identify  nodal  points  having  the  particular 
values  of  the  stream  function,  ,  corresponding  to  the 
streamline.  However,  except  in  rare  circumstances  the 
streamline  will  not  intersect  an  element  boundary  exactly 
at  a  nodal  point.  Therefore,  linear  interpolation  is  used 
to  determine  the  locations  at  which  the  streamlines 
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Figure  3.4  Determination  of  streamlines  from  numerical 
fluid  mechanics  solution 

intersect  the  element  boundaries.  After  scanning  every 
element,  a  collection  of  coordinates  is  obtained  which 
corresponds  to  locations  where  the  streamline  crosses  the 
element  boundaries.  These  points  must  next  be  placed  in 
order,  which  is  accomplished  by  noting  that  if  a  streamline 
crosses  an  element  boundary,  it  must  cross  a  second  bound¬ 
ary  in  the  same  element.  Thus,  given  the  initial  point  of 
the  streamline ,  shown  as  PQ  in  Figure  3.4,  the  neighbor  to 
this  point,  P^  ,  can  be  determined  since  it  lies  on  a 
boundary  to  the  same  element.  But  P^  is  also  a  point  of 
intersection  for  a  second  element.  In  fact,  all  stream 
points,  with  the  exception  of  the  first  and  last,  mark 
the  intersection  of  the  streamline  and  the  boundary  of 
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two  elements.  Having  found  the  second  element  containing 
P1  ,  the  next  neighbor,  P2  ,  may  be  found  since  it  lies 
in  the  same  element  as  P1  .  This  procedure  continues  un¬ 
til  all  the  stream  points  are  arranged  in  order ,  forming 
the  streamline. 

In  the  exceptional  case  where  a  streamline  crosses 
an  element  boundary  at  a  corner  node ,  the  stream  point 
marks  the  intersection  of  the  streamline  and  more  than  two 
elements  (see  Figure  3.5).  However,  simply  disregarding 
the  elements  which  the  streamline  does  not  pass  through 
results  in  the  stream  point  Joeing  contained  by  just  two 
elements . 

Having  formed  the  strealines,  it  remains  to 
numerically  determine  the  orientations  of  fibers  at  se¬ 
lected  locations  along  each  streamline.  To  do  this,  the 
coordinate  system  whose  axes  remain  fixed  in  direction  and 
translate  along  the  fiber  center  (the  x^  ,  x2  ,  x3  axes 
in  Figure  2.1)  is  adopted,  and  Equation  (3.3)  can  be  used 
to  determine  the  fiber  orientation  relative  to  this  sys¬ 
tem.  The  advantage  of  working  with  (3.3)  is  that  differ¬ 
entiation  is  with  respect  to  a  single  variable,  time.  To 
illustrate  the  procedure  for  integrating  (3.3),  suppose 
the  initial  fiber  orientation  is  given  at  point  PQ  in 
Figure  3.4  and  the  fiber  orientation  at  location  P^  is 
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Figure  3.5  Streamline  crossing  corner  node 

desired.  The  time  required  for  a  fiber  to  travel  from  PQ 
to  P^  can  be  estimated  by  dividing  the  distance  from  PQ 
to  P^  by  the  velocity  at  PQ  .  This  sets  the  time  increment 
for  evaluation  of  the  fiber  orientation  in  (3.3).  The 
vorticity  and  rate  of  deformation  components  are  approxi¬ 
mated  as  constants  over  the  length  of  the  element  equal 
to  those  values  computed  at  PQ  .  With  these  approxima¬ 
tions,  a  numerical  differential  equation  solver  using  the 
Runge-Kutta-Verner  fifth  and  sixth  order  method  is  adopted 
to  solve  (3.3).  The  differential  equation  solver  is  a 
subroutine  available  from  International  Mathematical  and 
Statistical  Libraries.  The  fiber  orientation  at  subse¬ 
quent  stream  points  is  obtained  in  an  analogous  manner. 
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The  stream  points  provide  an  effective  sampling  region  in 
which  to  plot  fiber  orientation  provided  that  the  elements 
are  appropriately  sized. 

Under  the  above  procedure  for  determining  fiber 
orientation,  it  is  obvious  that  the  finer  the  mesh  used, 
the  more  accurate  the  solution  provided  roundoff  errors  do 
not  become  significant.  Questions  of  convergence  of  this 
scheme  were  addressed  by  Givler  [5}  who  showed  that  this 
method  yields  sixth  order  convergence. 

3.4  Orientation  Parameters  for  Planar  Fiber 

Distribution 

To  this  point,  the  methodology  for  obtaining  the 
orientations  of  individual  fibers  at  discrete  locations  in 
the  flow  field  has  been  developed.  The  number  of  locations 
containing  prescribed  fiber  orientations  depends  on  the 
fineness  of  the  mesh  and  the  number  of  streamlines.  If 
fibers  of  varying  initial  orientations  are  input  at  the 
beginning  of  each  streamline  then  an  orientation  distri¬ 
bution  will  result  at  locations  downstream.  It  is  well 
known  that  the  material  properties  of  a  fiber-reinforced 
component  are  strongly  dependent  on  the  degree,  as  well 
as,  the  direction  of  fiber  orientation.  Since  the  ulti¬ 
mate  interest  is  in  determining  material  properties,  it 
is  necessary  to  develop  parameters  relating  the  degree  of 
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orientation  to  the  material  properties. 

In  a  given  region,  the  distribution  of  fibers  may 
be  described  by  an  orientation  density  function,  n.  For  a 
planar  distribution  of  fibers  (see  Figure  3.6),  the  orien¬ 
tation  density  is  a  function  of  the  single  angle,  <j>^ ;  viz, 

n  =  n  ( 4>1 ) 

The  orientation  density  function  is  periodic  with 
period  tt  since  the  angle  <}>^  +  tt  depicts  the  same 
orientation  ^ ;  hence, 

n  ( 4>x )  =  n  ( 

It  is  essential  to  confine  to  lie  within  an  interval 
equal  to  ir  which  isolates  one  period  in  the  orientation 
density  function. 

For  a  symmetric  orientation  distribution  about 
the  mode,  <}>£  ,  such  as  the  distribution  in  Figure  3.7,  one 
intuitively  expects  local  orthotropic  material  properties 
with  the  mode  angle  defining  a  principal  material  direc¬ 
tion.  The  validity  of  the  symmetric  distribution  assump¬ 
tion  is  investigated  for  a  real  flow  situation  in  Appendix 
1,  where  it  is  determined  that  the  symmetry  assumption 
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Xl 


Figure  3.6  Planar  fiber  distribution 

is  valid.  With  this  symmetry  the  orientation  density 
function  is  normalized  such  that 

(J>£+tt/ 2 

nt^Jd^  =  1  (3.6) 

*1 

In  an  interval  centered  at  the  mode,  the  mean,  and 

the  mode  are  equivalent. 

McCullough  et  al.  [13]  have  developed  a  technique 
for  relating  the  microstructure  in  a  short-fiber  composite 
material  to  its  material  properties,  which  utilizes  the 
concept  of  the  "aggregate  model."  The  microstructure  is 
modelled  as  a  collection  of  subregions,  or  grains,  of 
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Figure  3.7  Symmetric  orientation  distribution  about  the 
mode 

microlaminates  of  aligned  fibers.  McCullough  assumed  a 
symmetric  distribution  about  the  mode  and  introduced  a  new 
angle  £  which  is  measured  from  the  mode  (see  Figure  3.8). 


Figure  3.8  Orientation  density  as  a  function  of  5 
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Obviously,  <J>^  and  5  are  simply  related  through 


5  - 


(3.7) 


Under  this  model,  the  material  properties  are  related  to 

the  dispersion  through  two  parameters ,  f  and  g  ,  defined 

P  P 


f  =  2<cos  £>  -  1 
p 


(3.8) 


1  4 

gD  =  -j=-[8<cos  £>  -  3] 


(3.9) 


where  <cosra^> =  n(c)  cosmc  d<;.  The  orientation  param¬ 


eters  are  constructed  such  that  the  degree  of  collimation 

can  readily  be  interpreted.  For  f  =  g  =  1,  the  fibers 

P  P 

are  perfectly  aligned ,  whereas  f  =  g  =  0  signifies  a 

p  p 

random  distribution.  Intermediate  values  of  fp  and  gp 
represent  partial  degrees  of  orientation.  The  orientation 
parameters  may  be  expressed  in  terms  of  <j>.  by  substituting 
(3.6)  into  (3.9)  and  (3.10)  : 

4>  ^  '1-'TT  /  2 

f  =  2(  n(tj)1  )cos2  (4>_— 4>°)d4>-,  "1  (3.10) 

P  Jao  1  1  1 


<j)?+Tr/2 

1  f  1 
I  8 

J<f>£ 


nf^lcos  (4)^  — 3 


(3.11) 
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In  order  to  obtain  the  orientation  parameters  from  a  dis¬ 
tribution  in  cf> ^  ,  the  mode,  cj>£  ,  must  be  determined. 

In  this  work,  an  orientation  density  function  is 
not  explicitly  developed;  rather,  the  orientations  of  a 
finite  number  of  fibers  are  determined,  as  explained  pre¬ 
viously.  Hence,  the  orientation  parameters  must  be  ex¬ 
pressed  in  terms  of  a  finite  number  of  fibers;  this  is 
accomplished  in  the  next  section.  The  subsequent  section 
is  concerned  with  the  evaluation  of  (j>£  . 

3.4.1  Determination  of  the  Orientation  Parameters 

from  the  Orientations  of  a  Finite  Number  of  Fibers 


To  be  compatible  with  the  numerical  fiber 
orientation  scheme  developed  earlier,  the  orientation 
parameters  must  be  expressed  in  terms  of  the  orientations 
of  a  finite  number  of  fibers. 


To  determine  the  appropriate  expression  for  f  , 

ir 

one  begins  by  dividing  the  integral  in  (3.10)  into  a 
finite  number  of  intervals  and  approximating  the  behavior 
over  each  interval  via 


<|>?+Tr/2 

1  2 

n(4>1)cos  ((j>1-(j)°)d(J)1 

*1 


2 

s  Z  n  [  (<j>- )  .  ]  cos 
i=l  1  x 


[(4>1)i-4>|]A<()1 


(3.12) 


31 


where  k  is  the  number  of  intervals  and  (4^)^  is  the  average 

J.U 

value  of  4^  in  the  i  interval.  Defining  m^  as  the  num- 
her  of  fibers  whose  oreintations  lie  in  the  i  interval , 
it  can  be  determined  from  the  normalization  condition 
(3.6)  that 
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which  is  the  desired  expression  relating  f  to  the  orien 
tations  of  a  finite  number  of  fibers. 


Following  an  analogous  procedure,  the  following 
expression  for  g^  may  be  determined: 


P 


8  N  4 

—  Z  cos  (4).  -<J)° ) -3 

l_  j=l  j 


(3.16) 


The  expressions  (3.15)  and  (3.16)  are  summarized  in 
Table  3.1  for  ease  in  reference. 

With  the  above  definitions  of  the  orientation 
parameters, a  random  fiber  distribution  may  be  modelled 
from  a  finite  number  of  fibers  by  adjusting  the  fiber 
orientations  such  that  f  =  g  =  0 . 

p  3p 

Inspection  of  (3.15)  and  (3.16)  reveals  that  the 
mode,  4) |  ,  must  be  determined  to  compute  the  orientation 
parameters.  The  next  section  develops  a  method  for  deter¬ 
mining  the  mode  in  terms  of  the  orientations  of  a  finite 
number  of  fibers. 
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Table  3,1 

Orientation  Parameters  Expressed  in  Terms  of  the 
Orientation  of  a  Finite  Number  of  Fibers 


fp  =  1.2  OOS2!^  -*;)  -  1 

1—1  1 


gp  =  5 


N 


§  £  COS4  (<!>,  “<{>? )  “3 

Ni=l  1i  1 


3.4.2  Computation  of  the  Mode 


For  a  continuous  distribution  function,  the 
computation  of  the  mode  angle  is  quite  straightforward. 
One  simply  applies  the  definition  of  the  mode: 


dn  (4^) 

dt})^ 


0 


Alternatively,  in  the  case  cf  symmetric  distributions 
about  the  mode,  the  mode  is  equivalent  to  the  mean  in  an 
interval  centered  about  the  mode;  viz. 


>1 


(j)°+ir/2 

$1  n ( <j>.|  ) dcf)-  =  <4>1> 
J  4)£-ir/2 


(3.17) 
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The  condition  (3.17)  is  easier  to  apply  when  dealing  with 
distributions  described  in  terms  of  the  orientations  of  a 
finite  number  of  fibers.  The  procedure  to  determine  <f>£ 
using  (3.17)  is  illustrated  for  continuous  distributions 
and  then  extended  to  include  distribution  in  terms  of  a 
discrete  number  of  fibers. 

Consider  the  continuous  distribution  in  Figure 
3.7.  In  order  to  determine  the  mode  via  (3.17),  the 
period  centered  at  the  mode  must  be  established.  But 
this  interval  is  unknown  until  the  mode  itself  is  deter¬ 
mined.  Hence,  an  iterative  procedure  is  necessary  where 
a  mean  is  calculated  based  upon  an  initial  guess  for  the 
interval .  A  new  interval  centered  around  the  calculated 
mean  is  established,  from  which  a  new  mean  is  calculated, 
leading  to  a  new  interval.  The  procedure  continues  until 
the  difference  between  successively  calculated  mean 
values  lies  within  some  tolerance  requirement.  To  illus¬ 
trate  the  iterative  scheme ,  the  mode  for  the  distribution 
in  Figure  3.8  is  determined  iteratively  in  Figure  3.9. 

These  ideas  are  easily  extended  to  cover 
distributions  described  by  the  orientations  of  a  finite 
number  of  fibers.  Consider,  for  example,  the  distribution 
depicted  by  the  four  fibers  in  Figure  3.10.  Within  the 
range  0  <_  <j>^  <tt,  the  orientations  of  the  four  fibers  are 


N<Pl'3 


(c) 


<b)  (d) 


Figure  3.9  Determination  of  the  mode  via  the  iterative 
procedure.  In  (a)  the  mean  angle,  <<j)1>1,  is 
calculated  in  the  interval  (0,ir).  A  new  inter¬ 
val  centered  about  <$1^  is  established  in  (b) 
from  which  a  new  mean,  <<j>i>2'  is  determined. 
This  procedure  continues  until  a  close  estimate 
for  the  mode  is  established  in  (d) . 
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Figure  3.10  Determination  of  the  mode  from  the 

orientations  of  a  finite  number  of  fibers 
via  the  iterative  procedure 
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0 . 05tt  ,  0 . 15ir ,  0 . 25it  and  0.95it.  Calculating  the  mean,  <J>1 
from  the  standard  definition 


where  N  is  the  total  number  of  fibers  and  <J>1  is  the 

+*  Vi  ^ 

orientation  of  the  j  fiber, one  obtains  the  mean  angle 

depicted  in  Figure  3.10(a).  Centering  a  new  interval 
about  the  computed  mean  angle ,  the  orientation  of  one 
of  the  fibers  changes  from  0.95^  to  -0.05n  and  a  new  mean 
results,  which  is  shown  in  Figure  3.10(b).  A  further 
iteration  does  not  change  the  mean  value. 

3.5  Numerical  Solutions 

In  this  section  the  numerical  method  is  used 
to  solve  several  example  problems.  The  fiber  orienta¬ 
tion  in  plane  Poiseuille  flow  is  determined  first, 
and  the  results  are  compared  with  the  analytical  solution 
developed  in  Section  3.2.  Thereafter,  several  examples 
of  pragmatic  interest  are  presented. 

3.5.1  Plane  Poiseuille  Flow 


Consider  again  the  wall-bounded,  steady. 
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Newtonian  flow  in  Figure  3.1.  The  analytical  solution  for 
the  orientation  of  fibers  initially  perpendicular  to  the 
wall  boundaries  has  been  determined  in  Section  3.2.  In 
the  present  example,  the  fiber  orientation  is  determined 
numerically. 

As  previously  stated,  the  first  step  in  the 
numerical  scheme  is  to  solve  the  flow  equations  via 
the  finite  element  method.  Figure  3.11  depicts  the  finite 
element  mesh  and  associated  boundary  conditions  incorpo¬ 
rated  to  solve  the  flow  problem.  Newtonian  constitutive 
behavior  is  assumed. 

The  fluid  mechanics  results  are  presented  in  the 
form  of  contour  plots.  Stream  function,  pressure,  and 
vorticity  contours  are  plotted  in  Figures  3.12  through 
3.14,  respectively.  It  is,  of  course,  possible  to  calcu¬ 
late  these  quantities  analytically  for  this  simple  flow. 
Comparing  the  contour  plots  with  the  analytical  solutions 
reveals  that  the  numerical  flow  solution  is 
accurate. 

With  the  discretized  solution  to  flow  equations 
at  hand,  the  fiber  orientation  may  be  determined.  In 
Table  3.2,  a  comparison  is  made  between  the  numerical  and 
analytical  solutions  for  the  orientation  of  a  fiber  ini¬ 
tially  perpendicular  to  the  wall  boundary  traversing  the 


Figure  3.11  Finite  element  mesh  and  boundary  conditions  for  solution  of  plane 
Poiseuille  flow 


Figure  3.12  Stream  function  contours  (streamlines)  in  plane  Poiseuille  flow 
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re  3.14  Shear  stress  contours  in  plane  Poiseuille  flow 
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streamline  which  traces  the  path  =  0.6.  It  is  seen 
that  the  numerical  results  compare  favorably  with  the 
analytical  predictions. 


Table  3 . 2 

Comparison  of  Analytical  and  Numerical  Solutions 
in  Plane  Poiseuille  Flow  at  x„=0 . 6  for  a 
Fiber  Initially  Perpendicular  to  the 
Wall  Boundary 

4>n 

1  analytical ,  1  numerical, 

x3  rad  rad 


in 

• 

o 

-0.753 

-0.785 

1.0 

00 

o 

• 

rH 

1 

-1.11 

1.5 

-1.23 

-1.25 

It  is  of  interest  to  determine  contours  of  the 
mode  and  orientation  parameters  for  this  flow.  However, 
one  must  first  determine  the  number  of  fibers  needed  to 
adequately  portray  the  orientation  distribution.  This 
question  is  addressed  in  Appendix  2  where  it  is  determined 
that  as  few  as  ten  fibers  give  good  estimates  for  both 
the  mean  fiber  angle  and  orientation  parameters. 

A  fiber  plot,  shown  in  Figure  3.15,  reveals 
information  about  the  direction  of  orientation,  as 
well  as,  the  degree  of  fiber  alignment.  It  is  readily 
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Fiber  plot  for  initially  random  fiber  orientation  in  plane 
Poiseuille  flow 
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observed  from  the  figure  that  a  high  degree  of  alignment 
is  obtained  in  the  wall  region  with  the  direction  of 
collimation  parallel  to  the  wall.  Near  the  centerline, 
the  orientation  is  nearly  random.  These  observations  are 
reinforced  in  the  mode  angle  and  orientation  parameter 
contour  plots  presented  in  Figures  3.16  through  3.18. 

Fiber  alignment  in  the  wall  region  in  short-fiber  molded 
components  is  a  well  known  phenomenon. 

It  should  be  noted  that,  in  the  mode  angle  contour 

plot,  contours  are  not  plotted  in  highly  dispersed  regions 

(i.e.,  f  <0.4)  since  the  fibers  do  not  have  much  direction- 
P 

ality  in  these  highly  dispersed  areas.  The  contours  in 
the  more  oriented  regions  are  easier  to  interpret  since 
the  "noise"  imparted  by  plotting  the  contours  in  the  dis¬ 
persed  regions  has  been  eliminated.  This  condition  is 
applied  to  all  future  mode  angle  plots. 

3.5.2  Flow  Around  a  Circular  Inclusion  in  a  Finite 
Width  Channel 

A  problem  of  more  pragmatic  interest  is  the 
determination  of  fiber  orientation  in  the  region  around 
a  molded  hole  in  a  short-fiber  injection  molded  component. 
The  flow  geometry  consists  of  a  finite  width  channel  with 
a  centered  circular  inclusion.  In  this  example,  the  ratio 
of  channel  width  to  inclusion  diameter  is  10:1.  The 
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finite  element  mesh  and  boundary  conditions  to  determine 
the  flow  solution  are  depicted  in  Figure  3.19.  Newtonian 
behavior  is  assumed. 

The  boundary  conditions  require  explanation. 

Along  the  solid  boundaries  (i.e.  ,  the  wall  and  inclusion 
boundaries) ,  the  normal  no  slip  condition  is  assumed.  On 
the  centerline,  adjacent  to  the  inclusion,  the  shear  stress 
and  the  x2  component  of  velocity  are  both  equal  to  zero. 

At  the  input  and  exit  boundaries  of  the  flow,  it  is  assumed 
that  the  fluid  is  far  enough  away  from  the  inclusion  that 
the  flow  is  unaffected  by  it;  hence  Poiseuille  flow  exists 
along  these  boundaries . 

The  resulting  stream  function,  pressure,  and  shear 
stress  contours  are  plotted  in  Figures  3.20  through  3.22, 
respectively.  Here,  an  analytical  solution  is  not  avail¬ 
able,  thus  one  must  be  guided  by  intuition  in  interpreting 
the  fluid  mechanics  solution.  Inspection  of  the  contour 
plots  (especially  the  stream  function  contours)  reveal 
that  the  flow  solution  is  reasonable. 

One  next  determines  the  fiber  orientation.  Guided 
by  results  in  Appendix  2,  ten  randomly  oriented  fibers 
are  input  at  the  beginning  of  each  streamline.  The  re¬ 
sulting  orientation  parameter  and  mode  angle  contours  are 
depicted  in  Figures  3.23  through  3.25.  Several  interesting 
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observations  can  be  made  from  these  plots : 

1)  As  before  in  the  Poiseuille  flow  example, 

a  distinct  boundary  layer  of  aligned  fibers 
exists  in  the  region  near  the  wall. 

2)  A  second  boundary  layer  of  collimated  fibers 
is  present  resulting  from  flow  around  the 
insert.  Downstream  from  the  insert  there  is 
no  mechanism  available  to  misalign  the  fibers; 
consequently,  the  boundary  layer  propagates 
downstream. 

3)  A  core  of  random  fibers  exists  between  the 
two  boundary  layers. 

It  has  been  mentioned  previously  that  the  presence 
of  fibers  leads  to  pseudoplastic  fluid  behavior,  which  may 
be  portrayed  by  a  power  law  constitutive  assumption. 
Figures  3.26  through  3.34  present  orientation  parameter 
and  mode  angle  contours  for  various  power  law  indices 
under  unit  flow  conditions.  The  primary  observation  to 
be  made  is  that  lowering  the  power  law  index  leads  to  a 
larger  core  of  random  fibers. 

3.5.3  Flow  Around  a  Circular  Inclusion  in  an  Infinite 

Width  Channel 

In  the  previous  example,  the  wall  boundary  played 
an  important  role  in  determining  the  fiber  orientation. 
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Figure  3.19  Finite  element  mesh  and  boundary  conditions  for  determining  Newtonian 
flow  around  a  circular  inclusion  in  a  finite  width  channel 
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finite  width  channel 


Figure  3.21  Pressure  contours  for  Newtonian  flow  around  a  circular  inclusion 
in  a  finite  width  channel 
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Figure  3.22  Shear  stress  contours  for  Newtonian  flow  around  a  circular  inclusion 
in  a  finite  width  channel 
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Figure  3.25  Mode  angle  contours  resulting  from  Newtonian  flow  around  a  circular 
inclusion  in  a  finite  width  channel.  is  measured  clockwise  with 

respect  to  the  perpendicular  to  the  flow  direction. 
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om  flow  of  a  power  law  fluid  { index=0 . 8 ) 
in  a  finite  width  channel 
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Figure  3.28  Mode  angle  contours  resulting  from  flow  of  a  power  law  fluid 

(index=0.8)  around  a  circular  inclusion  in  a  finite  width  channel 
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Figure  3.29  Contours  of  fp  resulting  from  flow  of  a  power  law  fluid  (index-0.5) 
around  a  circular  inclusion  in  a  finite  width  channel 
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Figure  3.31  Mode  angle  contours  resulting  from  flow  of  a  power  law  fluid 

(index=0.5)  around  a  circular  inclusion  in  a  finite  width  channel 


ode  angle  contours  resulting  from  flow  of  a  power  law  fluid 
index=0.2)  around  a  circular  inclusion  in  a  finite  width  channel 
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One  may  wish  to  investigate  the  fiber  orientation  in  a 
channel  of  infinite  width,  thus  isolating  the  inclusion  as 
the  sole  mechanism  for  orienting  the  fibers. 

The  finite  element  mesh  and  associated  boundary 
conditions  to  determine  the  flow  field  are  presented 
in  Figure  3.35.  Newtonian  behavior  is  assumed.  This  ex¬ 
ample  presupposes  a  constant  unit  far-field  velocity  pro¬ 
file.  The  left, right,  and  top  wall  boundaries  are  placed 
far  enough  away  from  the  insert  that  the  far-field  velocity 
profile  exists  at  these  outer  boundaries.  The  resulting 
stream  function,  pressure,  and  shear  stress  contours  are 
depicted  in  Figures  3.36  through  3.38. 

Having  ascertained  the  flow  solution,  the  fiber 
orientation  is  determined.  Orientation  parameter  and  mode 
angle  contours  are  presented  in  Figures  3.39  through  3.41. 
These  figures  show  clearly  the  expected  presence  of  a 
boundary  layer  of  aligned  fibers  resulting  from  flow  over 


the  inclusion. 


0  0 
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Figure  3.36  Stream  function  contours  (streamlines)  for  Newtonian  flow  around 
a  circular  inclusion  in  an  infinite  width  channel 
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Figure  3.38  Shear  stress  contours  for  Newtonian  flow  around  a  circular 
inclusion  in  an  infinite  width  channel 


Figure  3.39  Contours  of  fp  resulting  from  Newtonian  flow  around  a  circular 
inclusion  in  an  infinite  width  channel 
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Figuree  3.40  Contours  of  gp  resulting  from  Newtonian  flow  around  a  circular 
inclusion  in  an  infinite  width  channel 
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Figure  3.41  Mode  angle  contours  resulting  from  Newtonian  flow  around  a  circular 
inclusion  in  an  infinite  width  channel.  is  measured  clockwise 

with  respect  to  the  perpendicular  to  the  flow  direction. 


CHAPTER  4.  AXISYMMETRIC  FLOW 


Having  developed  and  demonstrated  the  numerical 
method  for  determining  fiber  orientation  in  plane  flow, 
attention  is  now  turned  to  the  determination  of  fiber 
orientation  in  axisymmetric  flow.  Defining  a  cylindrical 
coordinate  system  as  shown  in  Figure  4.1,  it  is  known  from 
the  assumption  of  axisymmetry  that  the  9  direction  veloc¬ 
ity  as  well  as  the  variation  of  the  r  and  z  direction 
velocities  in  the  9  direction  of  the  fluid  are  identically 
equal  to  zero.  Consequently,  the  axisymmetric  flow  case 
represents  a  two-dimensional  flow  regime,  and  the  method 
developed  to  predict  fiber  orientation  in  plane  flow  is 
applicable  to  axisymmetric  flow.  However,  the  fiber 
orientation  equations  have  to  be  developed  under  the  axi¬ 
symmetric  simplifications.  This  is  accomplished  in 
Section  4.1.  Furthermore,  in  this  chapter,  the  assumption 
of  planar  fibers  is  lifted,  necessitating  the  development 
of  a  new  set  of  orientation  parameters.  This  topic  will 
be  discussed  in  Section  4.3.  The  analytical  solution  for 
fiber  orientation  in  Poiseuille  flow  is  determined  in 
Section  4.2  and  is  compared  with  the  numerical  solution 
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Figure  4.1.  Cylindrical  coordinate  system 


in  Section  4.4.  One  final  example,  the  fiber  orientation 
in  a  simulated  disk  molding  operation,  provides  a  practical 
illustration  of  the  numerical  solution  technique. 

4.1  Fiber  Orientation  Equations  in  Axisymmetric  Flow 

To  take  full  advantage  of  the  simplifications 
induced  by  axisymmetry,  modifications  must  be  introduced 
into  the  fiber  orientation  equations.  One  can  make  use 
of  the  axisymmetric  simplifications  if  the  fiber  orienta¬ 
tion  is  described  relative  to  a  cylindrical  reference 
frame  as  opposed  to  the  Cartesian  reference  frame  in 
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-3  - 

the  angle  between  the  fiber  and  the  y  (or  z)  axis.  To 
avoid  confusion,  it  is  renamed  as  B.  A  new  angle  a  is 
defined  as  the  angle  between  the  r  axis  and  the  projection 
of  the  fiber  on  the  x  -x  plane  in  a  plane  containing  the 
fiber  and  the  z  axis.  Reference  to  Figure  4.1  reveals 
that 


a  =  <j>^-9 


(4.1) 


where  0  is  the  cylindrical  coordinate  angle.  For  the 

analysis  that  follows,  the  cylindrical  coordinates  are 
-1  -2  -3 

defined  as  y  ,  y  ,  and  y  where 


Referring  to  Figure  4.2,  it  is  seen  that  the 
cylindrical  coordinates  (y1)  are  related  to  the  Cartesian 
coordinates  (x1)  by 


-2  -  7?  -1  -2 
x  =  rcosG  =  y  cosy 

-3  -  .  -  -1.-2 

x  =  rsmfl  =  y  siny 
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From  the  above  relations ,  the  following  partial  derivatives 
are  easily  evaluated: 


-2 

cosy 


-y^siny^ 


(4.2) 


-1  -2 
=  y  cosy 


0 


The  partial  derivative  terms  above  are  needed  for  tensor 
tranformations  which  will  be  performed  later  in  the  anal¬ 
ysis  . 

To  obtain  the  necessary  equations  to  solve  the 
fiber  orientation  in  axisymmetric  flow,  it  is  necessary 
to  transform  all  the  vorticity  and  rate  of  deformation 
components  in  Jeffery's  equations  (2.1)  and  (2.2)  to 
their  respective  components  in  cylindrical  coordinates. 

To  do  this ,  one  must  determine  the  tensor  character  of 
the  vorticity  and  rate  of  deformation  components  and  then 
transform  these  quantities  according  to  the  usual  rules 
of  tensor  transformation.  Having  transformed  the  tensor 
quantities ,  they  are  then  converted  into  physical 
components. 


Dealing  first  with  vorticity  components,  it  is 


known  by  definition  that 


zi  "  2eijktokj 


where  e. is  called  the  permutation  symbol  and  defined  by 
x  jk 


eijk 


r  0  if  any  two  indices  are  equal 

1  if  i,j,k  are  distinct  and  in  cyclic 
order  (i.e.,  123  or  231  or  312) 

-1  if  i,j,k  are  distinct  but  not  in  cyclic 
order  (i.e.,  132  or  213  or  321) 


Thus , 


Zi  =  h^32  ~  “23) 


z2  =  *<w13  ‘  W31: 


(4.3) 


z3  =  ^(W21  "  “l2} 


where  w  is  the  second  order  vorticity  tensor  defined  by 


“ij  “  %<vi.j  '  vj.i’ 


From  this  definition  it  is  obvious  that  =  -to^  and, 
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hence,  Equations  (4.3)  simplify  to 


Z1  “  "w23 


z2  ~  a13 


(4.4) 


z3  -  -^12 


Now,  the  vorticity  terms  in  (2.1)  and  (2.2)  are  Cartesian 
components  and  thus  can  represent  covariant,  contravariant 
or  physical  components  since  no  distinction  is  made  in  the 
Cartesian  system.  Since  the  partial  derivatives  (4.2)  are 
known,  the  vorticity  terms  are  chosen  to  transform  as  con¬ 
travariant  components  via 


18  ,  6xa  Sx^  -ij  . 

(x)  =  — r  — ^  u)  J  (y) 

<$y  Sy-1 


(4.5) 


Using  (4.2),  (4.4),  and  (4.5)  one  can  compute  z^ (x) 


z1 (x)  =  z1 (x) 


■u23(x) 


i|!  i§!  sil (J) 

6y  Sy-5 
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It  is  desirable  to  switch  to  physical  components,  which  is 
easily  accomplished  as  follows: 

^(x)  =  -r/g11/g^w<12>  (y) 

=  -'~<12>  (y ) 


where  g  is  the  metric  tensor. 

*\* 

Here,  the  usual  bracketed  notation  for  physcial 
components  is  eliminated  and,  instead,  the  indices  are 
subscripted  and  given  in  terms  of  the  coordinate  labels. 
The  same  procedure  is  employed  to  determine  z2(x)  and 
i3(x)  as  a  function  of  co<ij>(y),  and  the  results  are  pre¬ 
sented  in  Table  4.1. 
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The  Cartesian  components  of  the  rate  of  deformation 
tensor  are  computed  in  terms  of  the  cylindrical  components 
following  the  same  procedure.  First,  one  notes  that  the 
rate  of  deformation  tensor  is  second  order  and  defined  by 


d .  .  =  h  (v .  .  +v  .  . ) 

13  i.D  3/i 
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Utilizing  (4.2)  and  (4.6),  the  following  expression  may 
be  determined  for  ^2  (x)  : 


Again  converting  to  physical  components  via 


d<ij> 


’11 


(no  sum  on  i  or  j) 


one  can  determine 

d22(x)  =  cos20  drr  -  sin  29  drQ  +  sin29  dQ6 


The  transformation  of  the  remaining  rate  of  deformation 
components  follows  an  identical  procedure  and  the  results 
are  summarized  in  Table  4.2. 
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The  component  transformations  contained  in  Tables 
4.1  and  4.2,  as  well  as,  the  angle  changes  can  be  substi¬ 
tuted  into  (2.1)  and  (2.2)  to  obtain 

A2L  =  -oorQ+(sin0  (jor2+cos0  o>qz)  cos  (a+0)  cot8+ (sin9  to0z 

-cos0  torz)sin(a+0)cot8+B[  (sin0  d0z;-cos0  drz)sin(a 
+8)cot8+(hsin20(drr-d0e)+cos20  drQ }cos2 (a+0) + (sin0  drz ' 
+cos0  d0z) cos (a+0) cotB+M (d00-drr) 

+2sin20  drQ}sin2 (a+0) ]  (4.7) 


Table  4 . 2 

Cartesian  Components  of  Rate  of  Deformation  Tensor 
Expressed  in  Terms  of  Cylindrical  Components 


d^  (x)  dzz 

2-  -  -  -  2-  - 

d22(x)  =  cos  0  drr-sin20  drQ+sin  0  d0Q 

d33(x)  =  sin20  drr+sin20  dr0+cos20  dQ0 

d12(x)  =  cos 0  drz-sin8  dQz 

di3(x)  =  sin0  drz+cos0  dQz 

d23(x)  =  ^sin20(drr-d00)+cos20  drQ 
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41  =  (sine  u  +cos8  u>D  . )  sin  (a+0)  +  (-sine  wfl 
ot  rz  tiz 

+COS0  w  ) cos (a+0) +B[ (cose  d  -sin0  dfl  ) cos (a+0) cos26 
+^{%sin20 (drr-de0)+cos20  drQ}sin2 (a+0)sin28 
+  (sin0  drz+cos0  dQz)  sin  (a+0)  cos28+%{  (drr.-d00)  cos20 
-2sin20  drQ }cos2 (a+0) sin28+^(drr+de0 ) sin28]  (4.8) 

Since  axisymmetry  is  assumed,  (4.7)  and  (4.8)  are  indepen¬ 
dent  of  0,  and  can  be  simplified  by  inserting  any  discrete 
value  of  0,  which  results  in  the  following: 

=  -to  +5  cosacotB-^  sinacotB+B  [-d  sinacotB 
6t  r0  0z  rz  rz 

+dr0ccs2a+d0zcosacot8+Js  (50Q-3rr)  sin2a]  (4.9) 


4i  =  (joq  sina+o)  cosa+B[d  cosacos28+% (d  -d00)  cos2asin26 
6t  0z  rz  rz  rr 

+|(5rr+500)sin2S] 


(4.10) 
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These  equations  can  be  further  simplified  by  noting  that 
several  vorticity  and  rate  of  deformation  components  are 
identically  equal  to  zero  due  to  axisymmetry.  For  example, 


=  0 

Likewise,  it  can  be  shown  that  a)0z  ,  drQ  ,  and  dQz  are 
also  equal  to  zero.  Table  4.3  contains  a  complete  descrip¬ 
tion  of  the  physical  components  of  oa  and  d  in  terms  of  the 

'V  V 

velocities.  Cancelling  out  all  of  the  terms  which  are 
identically  zero  in  equations  (4.9)  and  (4.10), 

■fr  =  sinacotB+B  [-d  sinacotB+^s  (dfl0-d  )  sin2a]  (4.11) 

ot  rz  rz  oo  rr 

4r  =  co  cosa+B[5  cosacos2B+3s  (d  -dflfi)  cos2asin2B 
ot  rz  rz  rr 

t|(drr-de0)sin26]  (4.12) 

These  are  the  resulting  equations  describing  fiber  orien¬ 
tation  relative  to  cylindrical  coordinate  reference  axes 
in  axisymmetric  flow. 


To  this  point,  a  coordinate  system  which 
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translates  with  the  fiber  center  and  whose  axes  remain 
fixed  in  direction  has  been  adopted.  As  before,  one  may 
wish  to  describe  fiber  orientation  relative  to  an  inertial 
reference  frame.  A  prudent  choice  of  an  inertial  refer¬ 
ence  frame  is  one  whose  axes  lie  in  the  directions  of  the 

_ .  12  3.. 

y1  axes.  Thus,  the  cylindrical  frame  y  ,  y  ,  y  is  incor¬ 
porated  where  each  y^"  axis  of  this  system  lies  in  the  same 

-i  -1  -2  -3 

direction  as  the  corresponding  y  axis  of  the  y  ,  y  ,  y 
system.  In  this  inertial  reference  frame,  (4.11)  and 
(4.12)  transform  to 

^  =  -corzsinacot3+B  [-drzsinacot6+J5  (d0Q-drr)  sin2a]  (4.13) 

=  to  cosa+B[d  cosacos2S+%(d  -dfifi)cos2asin28 
Dt  rz  rz  rr  ott 

+f(drr+d00)sin23]  (4*14> 


where 


D  (  )  =  £1  ) 

Dt  <5t 


H  ) 

<5r 


+  v 


H  ) 

z  <5z 


is  the  substantial  time  derivative.  The  time  derivative 

6  [■  1  vanishes  when  the  boundary  conditions  are  time  inde 
ot 

pendent . 
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With  the  above  fiber  orientation  equations  at 
hand,  it  is  now  possible  to  determine  solutions  for  fiber 
orientation  in  axisymmetric  flow. 


Table  4 . 3 


Physical  Components  of  the  Rate  of  Deformation 


and  Vorticity  Tensors  in  Axisymmetric  Flow 


UKij>  =  3$ 
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4.2  Analytical  Solution  for  Fiber  Orientation  in 

Poiseuille  Flow 

The  pressure  driven,  steady,  laminar, incompressible 
flow  of  a  Newtonian  fluid  through  a  long,  smooth  round  tube 
provides  a  flow  situation  from  which  an  analytical  solution 
for  fiber  orientation  may  be  determined.  The  geometry  for 
this  flow  is  depicted  in  Figure  4.3.  For  simplicity,  a  unit 
radius  and  unit  flow  are  considered.  The  velocity 


Figure  4.3  Schematic  of  Poiseuille  flow 


components  are  independent  of  0 ,  thus  the  problem  is 
axisymmetric . 

The  flow  solution  is  not  difficult  to  determine 
since  the  flow  is  unidirectional: 
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vr  =  VQ  =  0 

v  =  2  (1-r2) 


Substituting  into  the  fiber  orientation  equations 
(4.13)  and  (4.14),  one  obtains  (for  time  independent 
boundary  conditions) : 


6a 

Sz 


r  sinacotS  (b-1) 
1-r2 


(4.15) 


66 

6z 


cosa 

2 


( 1-B  cos  26) 


(4.16) 


For  large  aspect  ratios  which  exist  in  fibers ,  B  is  very 
nearly  equal  to  unity  (i.e.;  for  r  =50,  B=2499/2501) . 
Consequently,  one  may  conclude  from  (4.15)  that 


6a 

6z 


0 
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Hence, a  is  a  constant  and  (4.16)  becomes  an  ordinary  dif 
ferential  equation  in  6.  For  fibers  initially  perpendic 
ular  to  the  flow  streamline  (i.e.,  3=tt / 2  at  z=0)  ,  the 
solution  to  (4.16)  is 


cot  6  =  r  tan 


-2r  z  cos  a 
(1-r2) (rp+l/rp) 


(4.17) 


The  solution  for  in-plane  fibers  is  plotted  in  Figure  4.4 
for  a  fiber  aspect  ratio  of  50.  A  comparison  of  Figures 
3.2  and  4.4  reveals  that  the  in-plane  fibers  show  identi¬ 
cal  orientation  behavior  with  in-plane  fibers  in  plane 
Poiseuille  flow.  This  is  an  expected  result  since  in¬ 
plane  fibers  experience  identical  fluid  deformations  and 
rotation  in  both  flows. 


Equation  (4.17)  provides  a  known  analytical 
solution  for  fiber  orientation  in  an  axisymmetric  flow 
although  it  is  necessarily  restricted  to  a  particular 
initial  fiber  orientation.  Nevertheless,  it  can  be  used 
to  check  the  validity  of  the  results  predicted  by  the 
numerical  scheme . 

To  complete  the  development  of  the  fiber 
orientation  scheme,  orientation  parameters  must  be  deter¬ 
mined  for  non-planar  orientation  states .  This  is  accom¬ 
plished  in  the  next  section. 
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Figure  4.4  Fiber  orientation  in  Poiseuille  flow  for  in-plane  fibers  initially 
perpendicular  to  flow  streamlines 
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4 . 3  Orientation  Parameters  for  Non-Planar  Fibers 

In  dealing  with  axisymmetric  flows,  the  planar 
fiber  orientation  assumption  has  been  lifted.  Thus  the 
orientation  distribution  becomes  a  function  of  two  angles 
as  shown  in  Figure  4.5;  viz, 

n  =  n  (y ,  c) 

such  that 


tt  ir/2 

n(y,?)siny  dy  d?  =  1 

o  o 

The  angle  c,  is  the  planar  angle  defined  in  Section  3.4 
while  y  is  an  axial  angle.  Both  y  and  £  are  measured 
from  principal  material  axes.  The  distributions  are 
assumed  to  have  the  following  symmetries : 

n(y,C)  =  n(y+7T,c) 
n(y,C)  =  n(-y,C+-rr) 
n(y,c)  =  n(y,-c) 


Here,  the  first  two  symmetries  reflect  the  fact  that  the 
orientation  angles  are  not  unique  and  the  last  symmetry 


Figure  4.5  Angle  descriptors  for  orientation  density 
function 


states  that  the  distributions  in  t,  is  symmetric  about  the 
2  axis.  If  one  further  assumes  that  the  orientation  den¬ 
sity  function  is  separable,  i.e., 

n (y, £)  =  n1 (y)n2 (c) 


then  McGee  [14]  has  shown  that  the  following  parameters 
relating  fiber  orientation  to  the  material  properties 
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result: 

TT  /  2 

'  2 

f  =  2  n0  U)cos  5  d;  -  1 

P  J  2 
o 

tt/2 

n2 (c)cos 


f 


a 


-  tt/2 

r  2 

3  cos  y  n1 (y) sinY  dY 

j 

-  o 


g 


a 


-  tt/2 

5  cos  y  n,  (ylsmy  dY 

-  o 


Here,  f  and  g  are  the  same  orientation  parameters  de- 
P  P 

fined  in  Equations  (3.8)  and  (3.9),  while  f&  and  ga 
represent  new  orientation  parameters  which  occur  in  the 
absence  of  the  planar  distribution  assumption. 

The  angles  y  and  £  must  be  related  to  the 
orientation  angles  a  and  3  in  Figure  4.2.  McGee  has  deter¬ 
mined  that  the  axial  angle  y  can  be  measured  from  any  of 
the  principal  material  axes.  If  the  z  axis  is  assumed  to 
be  a  principal  axis,  then  y  simply  represents  3.  With  the 
z  axis  chosen  as  one  principal  axis,  the  direction  asso¬ 
ciates  with  the  modal  value  of  a,  a0,  in  the  r-0  plane 
provides  an  estimate  for  a  second  principal  axis.  The 
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procedure  for  calculating  <j>|  presented  in  Section  3.4, 
may  also  be  utilized  to  determine  a°.  The  principal 
axes  are  depicted  in  Figure  4.6.  The  angle  z,  is  simply 
related  to  a  by 

z,  =  a-a° 


principal 

direction 


principal 

direction 

r 


Figure  4 . 6  Principal  axes  from  which  orientation 
distribution  is  measured 
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Incorporating  the  above  assumptions,  the 
orientation  parameters  may  be  expressed  in  terms  of  a  and 
8  via 


It  remains  to  determine  the  orientation  parameters 
in  terms  of  a  finite  number  of  fibers.  For  f  and  ,  the 
derivation  follows  identically  that  presented  in  Section 
3.4  and  the  following  expressions  result: 

2  N  2 

f  =  rr  I  cos^  (a4  -a0)  -  1 
P  Ni=1 

ifa  N  4 

g  siiz  cos  (a.-a°)-3 
P  5  Ni=i  1 

where  is  the  a-angle  associated  with  the  i 
N  is  the  total  number  of  fibers. 


fiber  and 
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Following  the  same  procedure,  the  following  expression  for 
g_  may  be  derived  : 

u 


100 


g 


a 


5 

N 


“  4 

-  £  cos  B . 

i=l  1 


1 


/  4 


The  expressions  relating  the  orientation  parameters  to  the 
orientation  of  individual  fibers  are  summarized  in  Table 

4.4. 


Table  4.4 

Orientation  Parameters  for  Non-Planar 
Fiber  Distributions 


2  N  2 

f  =  rz  l  cos^  (a.  -a0)  -  1 
P  Ni=i  1 


%  =  5 


4  £  cos^ (a--a°)  -  3 
Ni=l  1 


-C  — 


3  N  2 

4  £  cos  B-  -  1 

Ni=l  1 


12 


^a  = 


5  N  4 
-  £  cos  Bi  -  1 
i=l 


/  4 
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4.3.1  Random  Fiber  Orientation 

A  random  fiber  orientation  state  occurs  when  the 
orientation  parameters  are  all  equal  to  zero.  Table  4.5 
lists  the  orientation  angles  for  ten  fibers  which  lead  to 
a  random  orientation  distribution.  These  orientation 
angles  are  used  in  the  subsequent  examples  to  simulate  ini 
tially  random  fiber  distributions. 


Table  4.5 


Orientation  Angles  for  Ten  Randomly  Oriented  Fibers 


3  ,  rad 

0.314 

-0.565 

-0.723 

0.880 

-0.942 

1.100 

1.225 

-1.257 

1.382 

-1.414 


a,  rad 

0.942 

2.199 

0.314 

1.571 

2.827 

0.628 

1.885 

2.513 

0.0 

1.257 
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4 . 4  Numerical  Solutions 

In  this  section  numerical  solutions  for  fiber 
orientation  in  axisymmetric  flows  are  presented.  The  first 
example  involves  the  determination  of  the  fiber  orientation 
in  Poiseuille  flow,  and  the  results  are  compared  with  the 
analytical  solution  developed  in  Section  4.2.  In  a  second 
example,  the  fiber  orientation  in  a  simulated  disk  molding 
operation  is  determined. 

4.4.1  Poiseuille  Flow 

To  check  the  numerical  scheme,  the  numerical 
solution  for  fiber  orientation  in  Poiseuille  flow  (see 
Figure  4.3)  is  determined.  The  finite  element  mesh  with 
associated  boundary  conditions  to  solve  for  the  flow  is 
presented  in  Figure  4.7.  The  stream  function,  pressure 
and  shear  stress  contours  for  the  fluid  mechanics  solution 
are  shown  in  Figures  4.8  through  4.10. 

Having  ascertained  the  flow  solution,  the  fiber 
orientation  is  determined.  A  comparison  of  analytically 
and  numerically  predicted  orientation  is  presented  in 
Table  4.6  along  the  ip  =  9/16  streamline.  This  streamline 
corresponds  to  the  path  r=0 . 5 .  The  agreement  between  the 
numerical  and  analytical  values  is  excellent. 


4.0 
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It  is  also  of  interest  to  plot  the  orientation 
parameter  contours  for  an  initially  random  orientation  dis¬ 
tribution.  Reference  to  Section  4.2  reveals  that  the  a 
orientations  remain  constant;  hence,  one  need  not  plot  con¬ 
tours  of  a°,  f  ,  and  g  .  Plots  of  f  and  g  are  presented 
P  P  a.  d 

in  Figures  4.11  and  4.12,  respectively.  From  these  plots, 
one  detects  the  presence  of  the  boundary  layer  of  aligned 
fibers  parallel  to  the  wall  boundary  in  the  region  near 
the  wall  boundary. 


Table  4 . 6 

Comparison  of  Numerical  and  Analytical  Solutions 
for  Fiber  Orientation  in  Poiseuille  Flow  along 
the  Streamline  ip  =  9/16 


X 

Initial 

6, 

rad 

Orientations 

a, 

rad 

Numerical 

6, 

rad 

Solutions 

a, 

rad 

Analytical 

B, 

rad 

Solutions 

a  , 

rad 

0.5 

1.571 

0.0 

2.16 

0.0 

2.16 

0.0 

1.0 

1.571 

0.0 

2.50 

0.0 

2.50 

0.0 

1.5 

1.571 

0.0 

2.69 

0.0 

2.68 

0.0 

0.5 

1.571 

0.785 

2.02 

0.785 

2.01 

0.785 

1.0 

1.571 

0.785 

2.33 

0.785 

2.33 

0.785 

1.5 

1.571 

0.785 

2.54 

0.785 

2.53 

0.785 

60 
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Figure  4.11  Contours  of  fa  in  Poiseuille  Flow.  y  is  measured  with  respect  to 
the  flow  direction. 


110 


4.4.2  Simulated  Disk  Molding  Operation 

A  final  example  involves  simulating  an  actual 
axisymmetric  disk  molding  operation.  Figure  4.13  depicts 
the  geometry  for  the  problem.  In  the  usual  manner ,  the 
flow  solution  is  determined  first.  The  finite  element 
mesh  and  associated  boundary  conditions  for  Newtonian  flow 
are  presented  in  Figure  4.14.  A  constant  velocity  profile 
exists  along  the  boundary  where  the  fluid  enters  the  domain 
simulating  the  effect  of  a  plunger  pushing  the  fluid.  The 
usual  no-slip  condition  is  employed  along  the  wall  bound¬ 
aries.  The  normal  velocity  and  shear  stress  vanish  along 
the  centerline.  The  boundary  condition  along  the  exit 
boundary  is  the  analytically  determined  solution  for  radial 
flow  satisfying  the  continuity  and  Navier-Stokes  equations. 
The  resulting  stream  function,  pressure,  and  shear  stress 
contours  are  presented  in  Figures  4.15  through  4.17, 
respectively. 

With  the  flow  solution  at  hand,  the  fiber 
orientation  may  be  determined.  To  gain  a  better  visual 
perception  of  the  type  of  orientation  that  occurs ,  ini¬ 
tially  in— plane  fibers  (oc=0)  are  treated  first.  Examina¬ 
tion  of  (4.11)  reveals  that  in-plane  fibers  remain  in-plane; 
hence,  fiber  plots  can  be  drawn.  Figures  4.18  through  4.20 
present  fiber  plots  along  three  streamlines:  one  near 
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both  wall  boundaries  in  the  disk  section  and  one  which 
runs  through  the  disk  center.  Several  conclusions  may  be 
drawn  from  these  in-plane  fiber  plots: 

(1)  Rapid  alignment  of  in-plane  fiber  occurs  in 
the  sprue  section  of  all  three  streamlines. 
This  is  due  to  the  very  high  normal  stresses 
which  exist  in  the  converging  flow. 

(2)  In-plane  fibers  preferentially  orient 
parallel  to  the  flow  streamlines  near  the 
wall  boundaries  in  the  disk. 

(3)  In-plane  fibers  orient  perpendicular  to  the 
flow  streamlines  in  the  center  of  the  disk. 

Next,  the  orientation  is  determined  for  initial 
non-planar  random  orientation  distributions.  Figures  4.21 
through  4.25  present  the  resulting  contours  of  f&.  From 
these  figures,  the  following  conclusions  may  be  reached: 

(1)  A  very  high  degree  of  alignment  in  the 
z  direction  (f  >.0.9)  occurs  throughout 

cl 

the  converging  sprue  section. 

(2)  In  the  wall  region  in  the  disk,  the  fibers 
align  perpendicular  to  the  z  direction 
(fa^-0.4) 

In  the  center  section  of  the  disk,  the 

distribution  is  moderately  aligned  in  the 

z  direction  (f  -0.6)  near  the  sprue,  becoming 
a 


(3) 
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widely  dispersed  (f  2*0).  near  the  exit  region 
of  the  flow. 

Contours  of  f  are  shown  in  Figures  4.26  through 
P 

4.29.  From  these  figures,  one  may  conclude  that,  except 
for  a  strip  in  the  central  region  of  the  disk,  the  a 
alignment  is  nearly  perfect.  To  ascertain  the  direction  of 
cc  alignment,  contours  of  ct°  are  plotted  in  Figures  4.30 
and  4.31.  These  figures  reveal  an  in-plane  alignment  of 
fibers  (ct=0)  in  both  the  sprue  section  and  wall  regions. 

In  the  central  region  of  the  disk  a  moderate  degree  of 
orientation  exists  (f  =0.6)  with  a  varying  mode  angle. 

ci 

4.4.3  Comparison  of  Predicted  and  Actual  Orientation 

From  micrographs,  Ellery  [4]  has  determined  the 
fiber  orientation  in  axisymmetric  disks  under  slow  fill 
rates  for  glass  reinforced  phenolics.  Ellery  found  a  high 
(j0gj^0g  of  fiber  orientation  parallel  to  the  flow  stream— 
linos  in  the  wall  region  of  the  disk  and  a  high  degree  of 
orientation  in  the  "hoop"  direction  in  the  center  section 
of  the  disk.  This  experimentally  determined  orientation 
is  depicted  in  Figure  4.32. 

A  comparison  of  the  experimentally  observed 
orientation  with  the  numerical  prediction  reveals  a  very 
good  agreement  in  the  wall  region  of  the  disk.  In  the 
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center  section,  however,  the  numerical  scheme  does  not 
predict  a  high  degree  of  hoop  direction  reinforcement. 

Table  4.7  presents  the  orientation  of  ten  initially  ran¬ 
domly  oriented  fibers  at  the  beginning  of  the  disk  section 
and  at  the  exit  location  along  the  streamline  which  tra¬ 
verses  the  center  section  of  the  disk.  It  is  seen  that  the 
fiber  orientation  is  predominantly  in-plane  at  the  begin¬ 
ning  of  the  radial  section,  and  thus  hoop  direction  align¬ 
ment  cannot  occur.  However,  those  fibers  which  do  have  a 
significant  degree  of  out-of-plane  tilting  at  the  entrance 
to  the  disk  do  tend  to  orient  in  the  hoop  direction.  If 
one  introduces  a  fiber  tilted  45°  out-of-plane  at  the  en¬ 
trance  to  the  disk,  then  upon  reaching  the  exit  section 
the  fiber  has  achieved  a  hoop  direction  alignment.  Thus 
the  discrepancy  between  experimentally  observed  and  numer¬ 
ically  predicted  orentation  may  be  explained  by  postulating 
that  fiber  interactions  cause  a  significant  degree  of  out- 
of-plane  tilting  in  the  entrance  region  of  the  disk,  and 
the  orientation  mechanisms  are  able  to  orient  the  resulting 
non-planar  fibers  in  the  hoop  direction. 


Table  4.7 


Numerically  Predicted  Orientations  of  Individual 
Fibers  along  Center  Streamline  in  an 
Axisymmetric  Disk 


(3) 


(1) 

Initial 

Orientation, 

rad 

(2) 

Orientation  at  Begin¬ 
ning  of  Radial  Section, 
rad 

(3) 

Final 

Orientation 

rad 

3 

a 

6 

a 

3 

a 

0.314 

0.942 

-1.99 

3.1 

-2.79 

2.5 

-0.565 

2.199 

-1.92 

3.14 

-2.87 

3.10 

-0.723 

0.314 

2.05 

0.18 

2.38 

1.26 

0.880 

1.571 

-2.16 

3.13 

-2.87 

3.03 

-0.942 

2.827 

-2.50 

3.04 

-2.78 

2.34 

1.100 

0.628 

2.08 

0.19 

2.38 

1.27 

1.225 

1.885 

1.33 

3.12 

0.324 

2.75 

-1.257 

2.513 

2.18 

0.01 

2.87 

0.04 

1.382 

0.000 

2.18 

0.000 

2.87 

0.000 

-1.414 

1.257 

-1.82 

3.02 

-2.37 

1.89 

_ 

1.571 

0.785 

1.63 

1.54 

Figure  4.13  Schematic  for  molding  an  axisymmetric  disk 


116 


o 

0 

25 

256 

II 

ii  >i 

> 

Figure  4.14  Finite  element  mesh  for  axisymmetric  disk  molding  simulation 
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Figure  4.15  Stream  function  contours  (streamlines)  for  Newtonian  flow  in 
axisymmetric  disk  molding  simulation 
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molding  simulation 


Figure  '4.17  Shear  stress  contours  for  Newtonian  flow  in  axisymmetric  disk 
molding  simulation 
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Figure  4.18  In-plane  fiber  orientation  along  a  streamline  near  inside  wall  boundary 
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tr> 


re  4*19  In-plane  fiber  orientation  along  a  streamline  in  center  of  disk 
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wall 


23 


.21  Contours  of  f  =0.9  in  axisymmetric  disk  molding  simulation 
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Figure  4.22  Contours  of  f  =0.6  in  axisymmetric  disk  molding  simulation 
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Figure  4.23  Contours  of  f  =0.3  in  axisymmetric  disk  molding  simulation 
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.24  Contours  of  f  =0  in  axisymmetric  disk  molding  simulation 


Figure  4.25  Contours  of  f  =-0.4  in  axisymmetric  disk  molding  simulation 


Figure  4.27  Contours  of  f  =0.9  in  axisymmetric  disk  molding  simulation 
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Figure  4.28  Contours  of  f  =0.6  in  axisymmetric  disk  molding  simulation 


Figure  4.30  Contours  of  a°=0  in  axisymraetric  disk  molding  simulation 
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Figure  4.31  Contours  of  a°=1.0,  2.0,  3.0  in  axisymmetric  disk  molding  simulation 


CHAPTER  5.  CONCLUSIONS 


A  numerical  method  has  been  developed  for 
determining  the  fiber  orientation  in  arbitrary  plane  and 
axisymmetric  flows.  The  resulting  fiber  dispersion  is 
depicted  by  orientation  parameters ,  which  relate  the  degree 
of  fiber  alignment  to  the  material  properties.  In  axisym¬ 
metric  flow,  non-planar  distributions  are  treated  while 
only  planar  distributions  are  allowed  in  plane  flow.  The 
validity  of  the  numerical  method  was  proved  by  comparing 
numerical  results  with  known  analytical  solutions . 

Several  examples  were  considered.  In  Poiseuille 
flow,  a  distinct  boundary  layer  of  aligned  fibers  parallel 
to  the  flow  streamlines  was  shown  to  exist  in  the  wall  re¬ 
gion.  This  prediction  correlates  well  with  experimental 
observations .  It  was  determined  that  the  presence  of  a 
circular  inclusion  lead  to  boundary  layer  of  aligned  fibers 
adjacent  to  and  downstream  from  the  insert.  A  prediction 
for  the  fiber  orientation  in  an  axisymmetric  disk  was 
established.  The  prediction  correlated  well  with  experi¬ 
mental  observations  in  the  wall  region  of  the  disk  and  the 
discrepancy  in  the  central  region  could  be  explained  by 
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the  failure  of  the  theory  to  account  for  fiber  interactions. 

Several  extensions  to  the  theory  are  possible.  A 
more  elaborate  procedure  to  determine  the  principal  mate¬ 
rial  axes  in  non-planar  distributions  should  be  developed. 
In  the  axisymmetric  flow  examples  presented  in  Section  4.4, 
it  was  assumed  that  the  z  axis  represented  a  local  princi¬ 
pal  material  direction.  Clearly,  this  assumption  in  in¬ 
valid  for  many  flows,  and  a  procedure  needs  to  be  developed 
to  determine  the  principal  axes  from  the  orientations  of 
a  finite  number  of  fibers.  One  possible  method  of  attack 
is  to  assume  that  the  peak  in  the  orientation  density  func¬ 
tion  defines  one  principal  axis  and .  develop  a  pro¬ 
cedure  to  determine  the  peak.  Having  found  one  principal 
axis,  the  problem  is  no  more  difficult  than  the  determin¬ 
ation  of  the  principal  axes  in  a  planar  distribution. 

A  second  extension  involves  the  determination  of 
the  effect  of  the  fiber  orientation  on  the  viscosity  of 
the  suspension.  If  this  effect  is  significant,  then  an 
iterative  scheme  would  be  necessary  to  determine  the 
fiber  orientation,  as  duscussed  in  Section  2.1. 

Givler  [5]  has  discussed  several  other  possible 
extensions  to  the  theory,  including  studying  three- 
dimensional  and  transient  flows  and  determining  the 
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orientation  equations  for  a  viscoelastic  suspending  medium. 
These  extensions  enable  more  accurate  simulation  of  the 
molding  processes. 
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APPENDIX  1 


INVESTIGATION  OF  THE  SYMMETRY  OF  AN 
ORIENTATION  DISTRIBUTION 

In  developing  the  orientation  parameters,  the 
orientation  distribution  has  been  assumed  to  be  symmetric 
about  the  mode.  In  this  appendix,  the  validity  of  this 
assumption  is  investigated  for  planar  distributions  in 
the  plane  Poiseuille  flow  example  in  Section  3.5.1.  A 
random  fiber  orientation  distribution  is  input  at  the  be¬ 
ginning  of  the  \p=-0.29  streamline  and  the  skewness  of  the 
distribution  is  studied  at  selected  locations  downstream. 
The  degree  of  skewness  is  described  by  the  coefficient  of 
skewness  (a3)  defined  by 

a3  =  M3/(M2)3/2 

where  M2  and  M3  are  the  second  and  third  central  moments, 
respectively,  defined  by 

1  N 

Mi  =  n  £  (h-~  V1 
j=l  3 
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The  coefficient  of  skewness  is  determined  for  varying 
numbers  of  fibers. 

Figures  Al.l,  Al.2,  and  A1.3  depict  the 
orientations  of  100  initially  random  fibers  at  0.5,  1.0, 
and  1.5  inches  downstream,  respectively,  from  the  random 
end.  Also  presented  in  these  figures  are  the  numerically 
computed  mode  angles.  In  Figure  Al.4,  the  magnitude  of 
a^  is  plotted  as  a  function  of  the  number  of  fibers  at 
these  same  locations.  It  is  readily  observed  that  a^ 
decreases  to  small  values  as  the  number  of  fibers  is  in¬ 
creased.  Since  the  orientation  distribution  function  is 
more  accurately  modelled  by  larger  numbers  of  fibers,  it 
may  be  concluded  that  the  assumption  of  symmetric  distri¬ 
butions  is  valid  for  this  example.  While  the  conclusions 
from  this  specific  example  cannot  be  extended  to  other 
flows,  the  results  do  lend  credence  to  the  assumption 
of  symmetric  distributions. 
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vs.  number  of  fibers 


APPENDIX  2 


INVESTIGATION  OF  NUMBER  OF  FIBERS  NEEDED  TO 
ACCURATELY  PREDICT  ORIENTATION  PARAMETERS 

In  using  the  numerical  orientation  model,  one 
needs  an  estimate  of  how  many  fibers  are  needed  to  ade¬ 
quately  model  the  distribution  and  thus  give  accurate 
estimates  for  both  the  mode  angle  and  orientation  param¬ 
eters.  Obviously,  if  too  few  fibers  are  used,  grossly 
inaccurate  results  may  be  obtained.  On  the  other  hand, 
too  many  fibers  leads  to  wasted  computer  time. 

For  this  investigation,  orientations  were  studied 
at  selected  locations  along  the  ip=- 0.29  streamline  in  the 
plane  Poiseuille  flow  example  presented  in  Section  3.5.1. 
These  locations  are  identical  to  those  studied  in  Appendix 
1.  Plots  of  4>£  ,  fp  ,  and  g^  as  a  function  of  the  number 
of  fibers  are  presented  in  Figures  A2.1  through  A2.3, 
respectively,  for  locations  0.5,  1.0,  and  1.5  inches 
downstream  from  the  initial  point  in  the  streamline,  where 
a  random  orientation  state  is  input.  It  can  be  seen  from 
the  figures  that  all  parameters  attain  fixed  values  at  all 
locations  when  the  number  of  fibers  exceeds  50.  To 
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determine  the  dependence  on  numbers  of  fiber  for  fewer  than 
50  fibers,  magnified  plots  are  presented  in  Figures  A2.4 
through  A2.6.  From  these  figures,  one  can  conclude  that 
in  all  cases,  as  few  as  10  fibers  provide  accurate  predic¬ 
tions  for  all  the  parameters.  This  result  is  very  encour¬ 
aging  in  light  of  the  fact  that  it  is  desirable  to  work 
with  as  few  fibers  as  possible. 
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Figure  A2. 
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Figure  A2.4  Magnified  plot  of  <J>°  vs.  number  of  fibers 
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Figure  A2.5  Magnified  plot  of  f  vs.  number  of  fibers 
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Figure  A2.6  Magnified  plot  of  g  vs.  number  of  fibers 


